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Event-chain algorithm for the Heisenberg model: Evidence for 2 : ~ 1 dynamic scaling 
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We apply the event-chain Monte Carlo algorithm to the three-dimensional ferromagnetic Heisen¬ 
berg model. The algorithm is rejection-free and also realizes an irreversible Markov chain that 
satisfies global balance. The autocorrelation functions of the magnetic susceptibility and the energy 
indicate a dynamical critical exponent z « 1 at the critical temperature, while that of the magne¬ 
tization does not measure the performance of the algorithm. This seems to be the first report that 
the event-chain Monte Carlo algorithm substantially reduces the dynamical critical exponent from 
the conventional value of z ~ 2. 


I. INTRODUCTION 


Ever since the advent of the local Metropolis algo¬ 
rithm (LMC) [H, Monte Carlo simulations of systems 
with many degrees of freedom have played an impor¬ 
tant role in statistical physics. Near phase transitions, 
LMC is severely hampered by dynamical arrest phenom¬ 
ena such as critical slowing down for second-order transi¬ 
tions, nucleation and coarsening at first-order transitions, 
and glassy behavior in disordered systems. A number of 
specialized algorithms then allow to speed up the sam¬ 
pling of configuration space, namely the Swendsen-Wang 
Q and the Wolff Q cluster algorithms, the multicanoni- 
cal method Q and the exchange Monte Carlo method [1] 
based on extended ensembles. 

The above algorithms respect detailed balance, a 
sufficient condition for the convergence towards the 
equilibrium Boltzmann distribution. Recently, algo¬ 
rithms breaking detailed balance but satisfying the neces¬ 
sary global-balance condition have been discussed il. 
Among them, the event-chain Monte Carlo (ECMC) al¬ 
gorithm 0 has proven useful in hard-sphere ©[HI and 
more general particle systems [H, [l^ , allowing to equi- 
librate larger systems than previously possible [III fl¬ 
it has also been applied to continuous spin systems 15j. 
ECMC uses a factorized Metropolis filter [l^ and relies 
on an additional “lifting” variable to augment configu¬ 
ration space [fl. It is rejection-free and realizes an ir¬ 
reversible Markov chain. So far, however, the speedup 
realized by ECMC with respect to LMC has always rep¬ 
resented a constant factor in the thermodynamic limit, 
although larger gains are theoretically possible Bifl- 

In this paper, we apply ECMC to the three- 
dimensional ferromagnetic Heisenberg model, defined by 
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the energy 

Ei{S,}) = -jJ2Si-S„ ( 1 ) 

{iJ) 


where J is the unit of the energy. Si is a three-component 
unit vector and the sum runs over all neighboring pairs 
of the N = sites of a simple cubic lattice of linear ex¬ 
tension L. In our simulations, we consider the critical in¬ 
verse temperature /3c = T/Tc = 0.6930 [ll]. To describe 
the dynamics of the system, we compute the autocorre¬ 
lation functions of the energy, the system magnetization 
M = magnetic susceptibility 


_ |M|2 

“ N ■ 


( 2 ) 


The energy and the susceptibility are both invariant un¬ 
der global rotations of the spins Sk around a common 
axis, whereas the magnetization follows the rotation. We 
will argue that the energy and the susceptibility are slow 
variables, that is that their slowest time constant de¬ 
scribes the correlation (mixing) time of the underlying 
Markov chain. Under this hypothesis, we will present ev¬ 
idence that the ECMC for the three-dimensional Heisen¬ 
berg model reduces the dynamical critical exponent from 
the LMC value of z ~ 2 to z ~ 1. This considerable 
reduction of mixing times with respect to the LMC may 
well be optimal within the lifting approach [fl. The 
observed reduction is all the more surprising as in the 
closely related XY model [fl, where the spins are two- 
dimensional unit vectors, the ECMC realizes speedups 
by two orders of magnitude with respect to LMC, but 
does not seem to lower the dynamical critical exponent. 


II. ECMC ALGORITHM FOR THE 
HEISENBERG MODEL 

Applied to the Heisenberg model, the ECMC augments 
the physical space of spin configurations by a lifting vari¬ 
able (fc, v) which specifies the considered infinitesimal 
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FIG. 1. (Color online) Temperature dependence of the energy density e = E/N, the specific heat c and the magnetic suscepti¬ 
bility X of the three-dimensional Heisenberg model with L = 12. A chain length i = A^Tr/lO is used. 


counterclockwise rotation of spin k about the axis v. 
By virtue of the factorized Metropolis filter, this phys¬ 
ical move can only be rejected by a single neighboring 
spin, I, and the lifting variable will then be moved as 
{k,v) —>• {l,v), keeping the sense of rotation, but pass¬ 
ing it on to the spin responsible for the rejection. In 
the augmented space, the rejections are thus supplanted 
by events, namely the lifting moves for arrested physical 
states. The ECMC, for a given axis v, breaks detailed 
balance, yet satisfies global balance, as the probability 
flow into each lifted configuration equals the flow out of 
it, to first order in the time increment dt. For the XY 
model of planar rotators, v is uniquely defined as the 
axis perpendicular to the sense of rotation. For this rea¬ 
son, the FCMC around this axis is irreducible, and the 
chain length i in this model is best taken equal to the 
simulation time [lB |. For the Heisenberg model, spin ro¬ 
tations must be about at least two axes, in order to reach 
the entire configuration space. The resampling of the ro¬ 
tation axis is performed after the cumulative rotation 
angles about the previous axis reaches the chain length 
i. All configurations of the chain sample the equilibrium 
distribution and any uniform subset of them yield valid 
observable averages. Observables may be integrated dur¬ 
ing the continuous evolution or e.g. retrieved at regular 
intervals independent of the lifting events. 

For a fixed rotation axis the FCMC algorithm for 
the Heisenberg model reduces to the one of the XY 
model: With {(j)v,k,Sv,k) the spherical coordinates of a 
spin k in a system where the 2 -axis is aligned with v, the 
pair energy E^i between spins k and I is 

Eki — J cos(<)>ij^fc “t“ E (3) 


with 


angle dependence in Eq. m is oc cos{(pv,k — as in 

the XY model. 

The azimuthal angle (l)v,k increases for each ECMC 
chain from its initial value (j)o until one of its neighbors, 
I, triggers a lifting {k,v) -)■ {l,v) at = <()/, event- The 
latter is sampled with a sing le random number in the 
event-driven approach [i2j[l3]- Precisely, 0/,event is given 
by the sampling of the positive pair energy increase: 


AEi = - [log ran(0, 1)]/^ = 

j, dcofi{(t>^^k- (t>-c,i)\ 

-J max 0,-—- d0, 

J<t>0 V / 


(4) 


where ran(0,1) is a uniform random number between 0 
and 1. To solve Eq. (HD for 0/,event, One first slices off any 
full rotations (these n rotations by 2tt yield an energy 
increase of 2nJ'), leaving a value AF/, 

-^init + = - J' C0S((/);,event “ 'Pv,! “ SuTt), (5) 


where 



if initial pair energy change > 0 
otherwise. 


The true lifting event corresponds to the earliest of the 
independent event times sampled for all the neighbors of 
the spin k. In ECMC, Monte Carlo time is continuous 
and proportional to the total displacement of the spins. 
We have checked the correctness of the ECMC, and ob¬ 
tained perfect agreement for the mean energy, the specific 
heat and the susceptibility with the heat-bath algorithm 
dill modified with the exchange Monte Carlo method 
(or “parallel tempering”) @ (see Fig. [T|). 


J' = J sin^i, fc sin0„^;, 
K = — J cos9„^k cos 


III. DYNAMICAL SCALING EXPONENT 


Both J' and K depend only on the polar angles and At the critical temperature Tc, the correlation length 
remain unchanged along the event chain. The azimuthal- ^ of a model undergoing a second-order phase transition 
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FIG. 2. Autocorrelation functions and time constants of ECMC for the three-dimensional Heisenberg model at its critical 
point P = 0.693. Left: Energy density autocorrelation function Ge for for system sizes 4®, 8®, ..., 64®. Center: Susceptibility 
autocorrelation function Cx for ECMC for the three-dimensional Heisenberg system sizes 4®, 8®, .. ., 64®. Right: Scaling of the 
autocorrelation time Tx (resp. Te) of the susceptibility x (resp. energy density e) with system size L for ECMC {blue circles) 
(resp. {red triangles)) and of the autocorrelation time of the susceptibility for EMC {yellow squares). Error bars are smaller 
than the markers size. Right Inset: Speedup for the susceptibility x hr comparison to LMC for system sizes 4®, 8®, ..., 64®. 



FIG. 3. (Color online) Autocorrelation function of magneti¬ 
zation CM{t) at the critical temperature for various system 
sizes. The inset shows the spin autocorrelation function of a 
trivial algorithm that only performs global rotations in spin 
space along the two axes. 


equals the system size L and the autocorrelation time 
of slow variables r diverges as t ^ L^, where z is the 
dynamical critical exponent. We measure time in sweeps: 
One ECMC sweep corresponds, in average, to N lifting 
events and one LMC sweep to N attempted moves. Time 
autocorrelation functions are defined by 


Co{t) 


{0{t' + t)0{t'))-{0{t'))^ 
{0^{t')) - {0{t'))^ ’ 


( 6 ) 


where the brackets (• • •) indicate the thermal average and 
t' is set sufficiently large for equilibration. The dynam¬ 
ical critical exponent of LMC for the three-dimensional 
Heisenberg model was estimated from the autocorrela¬ 
tion function of the magnetization M as z = 1.96(6)[2l|. 


The over-relaxation algorithm (2^, seems to give 
2 ; ~ 1.10 [2lj which was obtained from the autocorre¬ 
lation function of the magnetization, and the Wolff algo¬ 
rithm is believed to yield a value close to zero: z > 0, 
a value obtained from the susceptibility autocorrelation 
function 

To evaluate the correlation time and the dynamical 
critical exponent for the ECMC, one must pay attention 
to the irreversible nature of the underlying Markov chain. 
During one event chain, spins all rotate in the same sense, 
and the system undergoes global rotations with taking 
into account the thermal fluctuation. This results in fast 
oscillations of the magnetization M and a quick decay 
of its autocorrelation function that is insensitive to the 
system size (see Fig. El, and even to the temperature. 
However, this effect is also visible for a trivial algorithm, 
which simply performs global rotations (see the inset of 
Fig. [3]). The trivial algorithm satisfies global balance, 
but its correlation time is infinite, as it does not relax the 
energy. A similar effect appears in the ECMC for particle 
systems , that likewise is not characterized by the mean 
net displacement of particles. To characterize the speed 
of the ECMC, we consider the energy density and the 
susceptibility that we conjecture to be slow variables at 
the critical temperature. Both x and e are insensitive to 
global rotations and do not oscillate. 

As shown in Fig. [21 the autocorrelation functions both 
of the energy density and of the susceptibility are well 
approximated as a single exponential decay 

Cx{t) =exp{-t/T) (7) 

on essentially the same timescales. Furthermore, the 
finite-size behavior of the autocorrelation times indicates 
z ~ 1 dynamical scaling. This z value is significantly less 
than for the LMC and very similar to the one obtained 
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for over-relaxation methods, although the ^ ~ 0 value of 
the cluster algorithm is not reached. 


IV. DISCUSSION AND SUMMARY 

The earliest application of lifting [T^ , the motion of a 
particle on a one-dimensional V-site lattice with periodic 
boundary conditions, already featured the decrease of the 
dynamical scaling exponent from z = 2 to 2 ; = 1 (the re¬ 
duction of the mixing time from oc N'^ to oc iV). To reach 
such reductions, the Markov chain must be irreversible. 
It was pointed out that the “square-root” decrease of the 
critical exponent was the optimal improvement [ 13 . The 
concepts of factorized Metropolis filters and of infinitesi¬ 
mal moves bring irreversible lifting algorithms to general 
fV-body systems, although only finite speed-ups were re¬ 
alized so far in the iV —^ 00 limit. The three-dimensional 
Heisenberg model seems to be a first such ECMC appli¬ 
cation with a lowered critical dynamical exponent. Our 
observation relies on the hypotheses that the energy and 
the susceptibility are indeed “slow” variables, and that 
the observed decay of the autocorrelation function con¬ 
tinues for larger times. However, in Fig. [21 a crossover 
from z = 1 back to z = 2 as it was observed in the XY- 
model after ~ 5 sweeps [l^ appears unlikely to arise af¬ 
ter hundreds of sweeps. The dynamical critical exponent 
z « 1 represents a maximal improvement with respect to 
the z « 2 of LMC, supposing again that the theorems of 
ref. [13 apply to infinitesimal Markov chains. 

In summary, we have successfully applied ECMC to 
the Heisenberg model in three dimensions. ECMC shows 


considerable promise for spin models, and the numerical 
data presented in this paper allow us to formulate the 
exciting conjecture that the dynamical critical exponent 
for the Heisenberg model is z ~ I. ECMC is also applica¬ 
ble to frustrated magnets and spin glasses, which involve 
antiferromagnetic interactions and/or quenched disorder. 
Our preliminary study indicates that the ECMC algo¬ 
rithm is also useful for a Heisenberg spin glass model. 
ECMC can be easily combined with other algorithms 
such as the exchange Monte Carlo method and the over¬ 
relaxation algorithm in the usual manner. This may al¬ 
low to investigate the three-dimensional Heisenberg spin 
glass model in the low-temperature region. Large-scale 
simulations in this direction are currently in progress. It 
would be very interesting to understand why ECMC is 
so much more successful in the Heisenberg model than 
both in hard and soft disks and in the XY model. 
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